\documentclass[letterpaper,10pt]{article}
\usepackage{geometry}   % See geometry.pdf to learn the layout
                        % options.  There are lots.
\usepackage[latin1]{inputenc}
\usepackage{graphicx}
\usepackage{epstopdf}
\usepackage{amsmath,amsfonts,amssymb}

\DeclareGraphicsRule{.tif}{png}{.png}{`convert #1 `dirname #1`/`basename #1 .tif`.png}

\author{Daniel R. Reynolds}
\title{{\tt gFLDSplit}: \\
A FLD-based Radiation and Chemistry Solver for ENZO}

\renewcommand{\(}{\left(}
\renewcommand{\)}{\right)}
\newcommand{\vb}{{\bf v}_b}
\newcommand{\xvec}{{\bf x}}
\newcommand{\Omegabar}{\bar{\Omega}}
\newcommand{\rhob}{\rho_b}
\newcommand{\dt}{\Delta t}
\newcommand{\Eot}{E^{OT}}
\newcommand{\Ef}{E_f}
\newcommand{\sighat}{\hat{\sigma}}
\newcommand{\Fnu}{{\bf F}_{\nu}}
\newcommand{\Pnu}{\overline{\bf P}_{\nu}}
\newcommand{\R}{I\!\!R}
\newcommand{\Rthree}{\R^3}
\newcommand{\eh}{e_h}
\newcommand{\ec}{e_c}
\newcommand{\Edd}{\mathcal F}
\newcommand{\Eddnu}{\Edd_{\nu}}
\newcommand{\mn}{{\tt n}}
\newcommand{\mB}{\mathcal B}
\newcommand{\mC}{{\mathcal C}}
\newcommand{\mL}{{\mathcal L}}
\newcommand{\mD}{{\mathcal D}}
\newcommand{\mDnu}{\mD_{\nu}}
\newcommand{\mCnu}{\mC_{\nu}}
\newcommand{\mLnu}{{\mathcal L}_{\nu}}
\newcommand{\mCe}{\mC_e}
\newcommand{\mLe}{\mL_e}
\newcommand{\mCn}{\mC_{\mn}}
\newcommand{\mLn}{\mL_{\mn}}


\textheight 9truein
\textwidth 6.5truein
\addtolength{\oddsidemargin}{-0.25in}
\addtolength{\evensidemargin}{-0.25in}
\addtolength{\topmargin}{-0.5in}
\setlength{\parindent}{0em}
\setlength{\parskip}{2ex}


\begin{document}
\maketitle

\section{Introduction}
\label{sec:intro}

This document describes a new, highly scalable,
field-based radiation and chemistry solver for Enzo, 
{\tt gFLDSplit}. The target applications of this solver include the
transport of radiation using a flux-limited-diffusion-based solver,
distributed among a possibly large number of 
processors.  In this solver, the radiation field is assumed to be
either a monochromatic radiation energy at the ionization threshold of
Hydrogen I ($h\nu = 13.6$ eV), or an integrated radiation energy
density with an assumed radiation spectrum.  This radiation field may
be coupled to the surrounding matter using either an assumption of
local thermodynamic equilibrium (no chemical ionization), or a model
for Hydrogen ionization. 

The defining characteristic of this solver in comparison
with the {\tt gFLDProblem} solver is that in a given time step, this
solver separately evolves the radiation and the gas energy
correction + primordial chemistry evolution system, using an
operator-split formalism.  The goal in using an operator split solver
on this problem is to provide for the most robust and computationally
efficient solver possible for these interacting processes.

This guide will only highlight the solvers and equations available in
this module.  For further details on the equations, numerical methods,
and verification tests relevant to this module, we refer to the paper 
\cite{ReynoldsHayesPaschosNorman2009}, that describes the related
fully-implicit version of this solver.  Moreover, we start right in
with the relevant parameters that must be defined to use the module.
The equations that describe what these options mean follow in
subsequent sections, with references placed accordingly.


\subsection{Current limitations}
\label{sec:limitations}

This solver is still a work in progress, though will serve many user's
needs in its current form.  Specific limitations that are present in
the current solver, and that are under active development include:
\begin{itemize}
\item Adaptive mesh refinement -- This solver may currently only be
  run in unigrid simulations.  This is because in the implicit
  radiation equation we must solve the problem on the entire mesh
  hierarchy, and the linear solver interface changes (and becomes more
  complex) when moving to a hierarchical mesh structure.  We have
  finished an initial interface for linear solvers on hierarchical
  meshes, in the context of self-gravity solves.  We plan to extend
  this interface to the implicit radiation solve within the 
  {\tt gFLDSplit} module in the near future. 
\item Increased chemistry -- Due to the use of a single radiation
  field with an assumed spectrum, our initial solver development
  focused only on primordial chemistry.  This was based on the fact that
  more intricate chemical models would be best simulated with a radiation
  approximation that allows spectral variation in space, to allow for
  spectral hardening and I-front preheating.  However, we do allow
  coupling of the {\tt gFLDSplit} solver to Enzo's chemistry and
  heating solvers, allowing investigation of problems in which
  spectral variation may be less important.
\end{itemize}




\section{{\tt gFLDSplit} usage}
\label{sec:module_usage}

In order to use the split FLD radiation solver module, and to
allow optimal control over the solver methods used, there are a number
of parameters that may be supplied to Enzo.  We group these into two
categories, those associated with the general startup of the module
via the Enzo infrastructure, and those that may be supplied to the
{\tt gFLDSplit} module itself.  However, prior to embarking on a
description of these parameters, there are a few requirements for any
problem that wishes to use the {\tt gFLDSplit} module. 

Foremost, Enzo must be built using the two configuration options
{\tt PHOTON} (enabled through the call {\tt gmake photon-yes} in the
Enzo source directory), and {\tt HYPRE} (enabled using the call 
{\tt gmake hypre-yes}).  Moreover, the machine Makefile must specify
how to include and link with an available HYPRE library (version $\ge$
2.8.0b).  If a user must compile HYPRE themselves to obtain this
version, they should make note of the HYPRE configuration option
{\tt --with-no-global-partition}, which must be used for solver
scalability when using over $\sim\!1000$ processors, but which results
in slower executables on smaller-scale problems.



\subsection{Startup parameters}

In a user's main problem parameter file, the following parameters
must be set (their default values are in brackets):
\begin{itemize}
\item {\tt RadiativeTransferFLD} [0] -- this must be set to 2.  Other
  values will either disable the FLD solver, or will use one in a
  non-desired context.
\item {\tt ImplicitProblem} [0] -- this must be set to 3 to use this
  {\tt gFLDSplit} module (among the possible implicit radiation solver
  modules).
\item {\tt ProblemType} [0] -- as usual, this is problem-dependent.
  However, for FLD-based solvers, the value of ProblemType
  should be within the 400's.
\item {\tt RadHydroParamfile} [NULL] -- this should contain the filename
  (with path relative to this parameter file) that contains all
  module-specific solver parameters (discussed below).  While the 
  {\tt gFLDSplit} module parameters may be supplied in the main
  parameter file, that filename must still be specified here (though
  it is not recommended, since the {\tt ReadParameterFile.C} routine
  will complain about all of the `unknown' parameters that are read
  elsewhere).
\item {\tt RadiationFieldType} [0] -- this can be any value {\em except}
  10 or 11, since those use pre-existing background radiation
  approximations.
\item {\tt RadiativeTransferFLDCallOnLevel} [0] -- this should currently
  be set to 0.  Future releases plan to enable the implicit solver on
  a statically nested subgrid, but this does not work at present.
\item {\tt RadiativeTransfer} [0] -- this must be set to 0.  A nonzero
  value will instead use an explicit ray-tracing solver for the
  radiation transport. 
\item {\tt RadiativeTransferOpticallyThinH2} [1] -- this must be set
  to 0.  A nonzero value will attempt to use a $1/r^2$ Lyman-Werner
  radiation field, that is {\em ignored} by the {\tt gFLDSplit} module.
\item {\tt RadiativeCooling} [0] -- this may be either 0 or 1.  A
  nonzero value will enable Enzo's built-in explicit subcycled gas
  cooling modules instead of the coupled radiative-ionization provided
  by {\tt gFLDSplit}. 
\item {\tt MultiSpecies} [0] -- this may be set to either 0 or 1.  A
  nonzero value will enable Enzo's built-in chemistry modules instead
  of the coupled radiative-ionization provided by {\tt gFLDSplit}.
\end{itemize}

In addition, if a user wishes to set up a new {\tt ProblemType} that
uses the {\tt gFLDSplit} module, they must allocate a standard Enzo
baryon field having the {\tt FieldType} set to {\tt RadiationFreq0}.
It is this baryon field that will be evolved by the {\tt gFLDSplit}
module, and that a user may access to obtain information on the
grey radiation field. 

Furthermore, the {\tt gFLDSplit} module currently computes its own
emissivity field $\eta(\xvec)$, in the routine \\
{\tt gFLDSplit\_RadiationSource.src90}.  A user may edit this file
to add in an emissivity field corresponding to their own 
{\tt ProblemType}.  Alternatively, a user may separately fill in the
BaryonField {\tt Emissivity0} in some other Enzo routine.  Then, when
Enzo is compiled using the pre-processor directive {\tt EMISSIVITY}
and run with the global parameter {\tt StarMakerEmissivity $\ne$ 0},
the {\tt gFLDSplit} module will use the {\tt Emissivity0} field
instead of computing its own. 



\subsection{Module parameters}

Once a user has enabled the {\tt gFLDSplit} module, they have
complete control over a variety of internal module parameters.  The
parameters are given here, with their default values specified in
brackets, and references to the appropriate equations elsewhere in
this document. 
\begin{itemize}
\item {\tt RadHydroESpectrum} [1], this parameter chooses the type of
  assumed radiation energy spectrum from equation \eqref{eq:spectrum}.
  Allowed values include
  \begin{itemize}
  \item[1.] $T=10^5$ K blackbody spectrum, 
    \[
       \chi(\nu) = \frac{8 \pi h
         \left(\frac{\nu}{c}\right)^3}{\exp\left(\frac{h\nu}{k_b 10^5}\right)-1}.
    \]
  \item[0.] power law spectrum,
    \[
      \chi(\nu) = \left(\frac{\nu}{\nu_{HI}}\right)^{-1.5}
    \]
  \item[-1.] monochromatic spectrum at frequency $h\nu_{HI} = 13.6$ eV.
  \item[-2.] monochromatic spectrum at frequency $h\nu_{HeI} = 24.6$ eV.
  \item[-3.] monochromatic spectrum at frequency $h\nu_{HeII} = 54.4$ eV.
  \end{itemize}
\item {\tt RadHydroChemistry} [1], this parameter controls how {\tt
  gFLDSplit} interfaces with chemistry.  While Enzo may be run with
  any number of allowable chemical species, {\tt gFLDSplit} will
  interact with primordial chemistry through opacities, photo-heating
  and photo-ionization.  Allowable values are
  \begin{itemize}
  \item[0.] no chemical interactions
  \item[1.] Hydrogen chemistry
  \item[3.] Hydrogen + Helium chemistry
  \end{itemize}
\item {\tt RadHydroHFraction} [1], this parameter controls the
  fraction of baryonic matter comprised of Hydrogen, allowable
  values are $0 \le {\tt RadHydroHFraction} \le 1$.
\item {\tt RadHydroModel} [1], this parameter determines which model
  for radiation-matter coupling we wish to use, allowable values
  include
  \begin{itemize}
  \item[1.] Use the chemistry-dependent model from section
    \ref{sec:chem_model}, with a case B HII recombination coefficient.
  \item[4.] The same as model 1, but assume an isothermal gas energy
    (for regression test problems).
  \item[10.] Use the local thermodynamic equilibrium model from section
    \ref{sec:lte_model}.
  \end{itemize}
\item {\tt RadHydroMaxDt} [$10^{20}$], this parameter sets the value of
  $\dt_{\text{max}}$ from section \ref{sec:dt_selection}; it must be
  greater than 0.  This value must be given in {\em scaled} time units,   
  i.e.~$\dt_{\text{physical}} \le \dt_{\text{max}}*\text{TimeUnits}$, 
  where TimeUnits is Enzo's internal time scaling factor for the simulation.
\item {\tt RadHydroMinDt} [0], this parameter sets the value of
  $\dt_{\text{min}}$ from section \ref{sec:dt_selection}; it must be
  non-negative.  This value must also be given in scaled time
  units.
\item {\tt RadHydroInitDt} [$10^{20}$], this parameter sets the initial
  time step size for the {\tt gFLDSplit} module.  We note that since
  the module will take the smaller of $\dt_{\text{FLD}}$ and
  $\dt_{\text{CFL}}$, the default value is never actually used.  This
  value must also be given in scaled time units.
\item {\tt RadHydroMaxSubcycles} [$1$], this parameter sets the
  desired number of {\tt gFLDSplit} time steps per each hydrodynamics
  time step (must be $\ge 1$).  When using Enzo's chemistry and
  cooling solvers this parameter should be set to 1 to avoid overly
  decoupling radiation and chemistry.
\item {\tt RadHydroMaxChemSubcycles} [$1$], this parameter sets the
  desired number of chemistry and heating time steps per each
  radiation time step, {\em within} the {\tt gFLDSplit} solver.  This
  value is unused if Enzo's chemistry and heating routines are used.
\item {\tt RadHydroDtNorm} [2], this parameter sets the value
  of $p$ from equation \eqref{eq:time_error}.  
  \begin{itemize}
  \item A value of $0$ implies to use the $\max$ norm, 
  \item A value $>0$ implies to use the corresponding $p$-norm,
  \item Values $<0$ are not allowed.
  \end{itemize}
\item {\tt RadHydroDtGrowth}, [$1.1$], this value sets the maximum
  growth factor in the {\tt gFLDSplit} time step size from one time
  step to the next.
\item {\tt RadHydroDtRadFac}, {\tt RadHydroDtGasFac} and {\tt
    RadHydroDtChemFac}  [$10^{20}$], these parameters give the values
  of $\tau_{\text{tol}}$ for the variables $E$, $e_c$ and $\mn_{HI}$
  from equation \eqref{eq:time_estimate}, respectively.  They must be
  positive; the default specifies no restrictions on $\dt_{\text{FLD}}$.
\item {\tt RadiationScaling}, {\tt EnergyCorrectionScaling} and 
  {\tt ChemistryScaling} [1.0], these parameters give the scaling 
  factors $s_E$, $s_e$ and $s_{\mn}$ from
  \eqref{eq:variable_rescaling}, respectively; supplied values must be
  positive. 
\item {\tt AutomaticScaling} [1], this enables an heuristic approach
  to set the above scaling factors automatically by {\tt gFLDSplit}.
  This has been heavily utilized within reionization calculations.
  However, this is not recommended if optimal scaling factors are
  known for specific problems.
\item {\tt RadHydroTheta} [1.0], this parameter specifies the
  value of $\theta$ in equation \eqref{eq:radiation_PDE_theta},
  $0\le\theta\le 1$.
\item {\tt RadiationBoundaryX0Faces}, {\tt RadiationBoundaryX1Faces}
  and {\tt RadiationBoundaryX2Faces} [0 0], these specify the
  boundary-condition types from section \ref{sec:boundary_conditions}
  to use on the lower and upper boundaries in each direction.
  Allowable values are
  \begin{itemize}
  \item[0.] periodic (must match on both faces in a given direction)
  \item[1.] Dirichlet
  \item[2.] Neumann
  \end{itemize}
\item {\tt RadHydroSolTolerance} [$10^{-8}$], this parameter
  specifies the linear tolerance $\delta$ from section
  \ref{sec:rad_solve}.  Allowable values must be between $10^{-15}$ and 1.
\item {\tt RadHydroKrylovMethod} [$1$], this parameter
  specifies the type of outer linear solver algorithm.  Allowable
  values are
  \begin{itemize}
  \item[0.] Preconditioned Conjugate Gradient (PCG)
  \item[1.] Stabilized Bi-Conjugate Gradient (BiCGStab)
  \item[2.] Generalized Minimum Residual (GMRES)
  \end{itemize}
\item {\tt RadHydroMaxMGIters} [50], this positive parameter
  specifies the maximum number of multigrid iterations to perform in
  the MG-CG solver from section \ref{sec:rad_solve}.
\item {\tt RadHydroMGRelaxType} [1], this parameter specifies the
  relaxation method used by the multigrid solver:
  \begin{itemize}
  \item[0.] Jacobi
  \item[1.] Weighted Jacobi
  \item[2.] Red/Black Gauss-Seidel (symmetric)
  \item[3.] Red/Black Gauss-Seidel (nonsymmetric)
  \end{itemize}
  For more information, see the HYPRE user manual.
\item {\tt RadHydroMGPreRelax} [1], this positive parameter
  specifies the number of pre-relaxation sweeps the multigrid solver
  should use in the MG-CG solver from section \ref{sec:rad_solve}.
\item {\tt RadHydroMGPostRelax} [1], this positive parameter
  specifies the number of post-relaxation sweeps the multigrid solver
  should use in the MG-CG solver from section \ref{sec:rad_solve}.
\item {\tt EnergyOpacityC0}-{\tt EnergyOpacityC2} [1, 1, 0],
  these specify the opacity-defining constants $C_0$-$C_2$ for the
  energy-mean opacity $\kappa_E$ as in equation \eqref{eq:opacityE}.
  The parameters $C_0$ and $C_2$ must be $\ge 0$, while $C_1> 0$.
\end{itemize}





\section{Flux-limited diffusion radiation model}
\label{sec:rad_model}

We begin with the equation for flux-limited diffusion radiative
transfer in a cosmological medium \cite{ReynoldsHayesPaschosNorman2009},
\begin{equation}
\label{eq:radiation_PDE}
  \partial_{t} E + \frac1a \nabla\cdot\(E\vb\) =
    \nabla\cdot\(D\,\nabla E\) - \frac{\dot{a}}{a} E - c\kappa E + \eta,
\end{equation}
where here the comoving radiation energy density $E$, emissivity
$\eta$ and opacity $\kappa$ are functions of space and time.  In this
equation, the frequency-dependence of the radiation energy has been
integrated away, under the premise of an assumed radiation energy
spectrum, 
\begin{align}
  \notag
  & E_{\nu}(\nu,\xvec,t) = \tilde{E}(\xvec,t) \chi(\nu), \\
  \notag
  \Rightarrow & \\
  \label{eq:spectrum}
  & E(\xvec,t) = \int_{\nu_{HI}}^{\infty} E_{\nu}(\nu,\xvec,t)\,\mathrm{d}\nu 
    = \tilde{E}(\xvec,t) \int_{\nu_{HI}}^{\infty} \chi(\nu)\,\mathrm{d}\nu,
\end{align}
where $\tilde{E}$ is an intermediate quantity (for analysis) that is
never computed.  We note that if the assumed spectrum is the Dirac
delta function, $\chi(\nu) = \delta_{\nu_{HI}}(\nu)$, $E$ is a
monochromatic radiation energy density at the ionization threshold of
HI, and the $-\frac{\dot{a}}{a}E$ term (obtained through integration
by parts of the redshift term
$\frac{\dot{a}}{a}\partial_{\nu}E_{\nu}$) is omitted from
\eqref{eq:radiation_PDE}. Similarly, the emissivity function
$\eta(\xvec,t)$ relates to the true emissivity 
$\eta_{\nu}(\nu,\xvec,t)$ by the formula
\begin{equation}
\label{eq:emissivity}
  \eta(\xvec,t) = \int_{\nu_{HI}}^{\infty}\eta_{\nu}(\nu,\xvec,t)\,\mathrm{d}\nu.
\end{equation}

The function $D$ in the above equation \eqref{eq:radiation_PDE} is
the {\em flux-limiter} that depends on $E$, $\nabla E$ and the 
opacity $\kappa$,  
\[
   D(E) = \text{diag}\( D_1(E),\, D_2(E),\, D_3(E) \),
\]
where the directional limiters $D_i(E)$ are given by \cite{Morel2000}
\begin{align}
  \label{eq:Larsen_limiter}
   D_i(E) = \frac{c}{\sqrt{(3\kappa)^2 + R_i^2}}, \qquad
   R_i(E) = \max\left\{\frac{|\partial_i E|}{E}, 10^{-20} \right\}.
\end{align}




\section{LTE couplings}
\label{sec:lte_model}

For problems in which chemical ionization is unimportant, we may
assume that the gas is in local thermodynamic equilibrium.  In this
case we do not need to couple the radiation energy to a model for
chemical ionization, we therefore couple the radiation to the specific
gas energy equation, 
\begin{align}
  \label{eq:cons_energy}
  \partial_t e + \frac1a\vb\cdot\nabla e &=
    - \frac{2\dot{a}}{a}e
    - \frac{1}{a\rhob}\nabla\cdot\left(p\vb\right) 
    - \frac1a\vb\cdot\nabla\phi + G - \Lambda.
\end{align}
All but the final two terms in \eqref{eq:cons_energy} are already
handled by Enzo's existing hydrodynamics solver infrastructure.  In
this module, we therefore consider a specific energy correction equation 
\begin{align}
  \label{eq:cons_energy_correction}
  \partial_t e_c &= -\frac{2\dot{a}}{a}e_c + G - \Lambda,
\end{align}
that we use to correct Enzo's original gas energy to include radiation
couplings.  Here $G$ is the local heating rate,
\begin{align}
\label{eq:G_LTE}
  G &= \frac{c \kappa}{\rhob} E,
\end{align}
and $\Lambda$ corresponds to the local cooling rate,
\begin{align}
\label{eq:Lambda_LTE}
  \Lambda = \frac{\eta}{\rhob} .
\end{align}
The user-defined energy mean opacity $\kappa$ is
spatially-homogeneous, and is given by the formula
\begin{align}
\label{eq:opacityE}
  \kappa = C_0 \left(\frac{\rhob}{C_1}\right)^{C_2}
\end{align}
(the constants $C_0\to C_2$ are input by the user), and $\eta$ is
a black-body emissivity given by 
\begin{align}
\label{eq:etaBB}
  \eta = 4\kappa\,\sigma_{SB}\,T^4,
\end{align}
where $\sigma_{SB}$ is the Stefan-Boltzmann constant [$5.6704\times
10^{-5}$ erg s$^{-1}$ cm$^{-2}$ K$^{-4}$], and $T$ is the gas
temperature [K]. 



\section{Chemistry-dependent couplings (internal)}
\label{sec:chem_model}

In general, radiation calculations in Enzo are used in simulations
where chemical ionization states are important.  For these situations, 
we couple the radiation equation \eqref{eq:radiation_PDE} with
equations for both the specific gas energy correction and the
ionization dynamics of Hydrogen,
\begin{align}
  \notag
  \partial_t e_c &= -\frac{2\dot{a}}{a}e_c + G - \Lambda, \\
  \label{eq:hydrogen_ionization}
  \partial_t \mn_{HI} + \frac{1}{a}\nabla\cdot\(\mn_{HI}\vb\) &=
    \alpha^{rec} \mn_e \mn_{HII} - \mn_{HI} \Gamma_{HI}^{ph}. 
\end{align}
Here, $\mn_{HI}$ is the comoving Hydrogen I number density.  The
recombination rate $\alpha^{rec}$ is given by the case-B
recombination rate, 
\begin{equation}
\label{eq:alphaB}
\alpha^{rec} = 2.753\times 10^{-14} \left(\frac{3.15614\times 10^5}{T}\right)^{3/2} 
                   \left(1+\left(\frac{3.15614\times 10^5}{2.74\, T}\right)^{0.407}\right)^{-2.242}.
\end{equation}

In this model, the gas heating and cooling rates are
chemistry-dependent, 
\begin{align}
  \label{eq:G_nLTE}
  G &= \frac{c\,E\,\mn_{HI}}{\rhob} 
    \left[\int_{\nu_{HI}}^{\infty} \sigma_{HI}\, \chi_E
    \left(1-\frac{\nu_{HI}}{\nu}\right)\, d\nu\right] \bigg/
    \left[\int_{\nu_{HI}}^{\infty} \chi_E d\nu\right], \\
\label{eq:Lambda_nLTE}
  \Lambda &= \frac{\mn_e}{\rhob}\bigg[\text{ce}_{HI}\, \mn_{HI} 
  + \text{ci}_{HI}\, \mn_{HI} + \text{re}_{HII}\, \mn_{HII} + \text{brem}\,
  \mn_{HII} \\
  \notag &\qquad+ \frac{m_h}{\rho_{units}\, a^3} \left(\text{comp}_1\, (T-\text{comp}_2) 
    + \text{comp}_{X}\, (T-\text{comp}_{T})\right) \bigg].
\end{align}
The temperature-dependent cooling rates
$\text{ce}_{HI}$, $\text{ci}_{HI}$, $\text{re}_{HII}$, $\text{brem}$,
$\text{comp}_1$, $\text{comp}_2$, $\text{comp}_{X}$ and
$\text{comp}_{T}$ are all taken from Enzo's built-in rate tables.

The photo-ionization rate is calculated as
\begin{align}
  \label{eq:phHI}
  \Gamma_{HI}^{ph} \ = \ \int_{\nu_{HI}}^{\infty}
    \frac{c\,E_{\nu}\,\sigma_{HI}}{h\nu}\, d\nu 
  \ = \ \frac{c\,E}{h} 
    \left[\int_{\nu_{HI}}^{\infty} \frac{\sigma_{HI}\, \chi_E}{\nu}\,
      d\nu\right] \bigg/ \left[\int_{\nu_{HI}}^{\infty} \chi_E d\nu\right].
\end{align}
Similarly, the frequency-integrated opacity is now chemistry-dependent,
\begin{equation}
\label{eq:opacityHI}
  \kappa \ = \ 
  \left[\int_{\nu_{HI}}^{\infty} \kappa_{\nu}\,E_{\nu}\,d\nu\right] \bigg/
  \left[\int_{\nu_{HI}}^{\infty} E_{\nu}\,d\nu\right] \ = \ 
  \mn_{HI} \left[\int_{\nu_{HI}}^{\infty}
    \chi_E\,\sigma_{HI}\,d\nu\right] \bigg/
  \left[\int_{\nu_{HI}}^{\infty} \chi_E\,d\nu\right],
\end{equation}
where these integrals with the assumed radiation spectrum $\chi(\nu)$
handle the change from the original frequency-dependent radiation
equation to the integrated grey radiation equation.




\section{Chemistry-dependent couplings (external)}
\label{sec:chem_model_enzochem}

The primary operating mode for the {\tt gFLDSplit} solver is to couple
the FLD solver to Enzo's chemistry and heating approach.  To support
these situations, we couple the radiation equation
\eqref{eq:radiation_PDE} with Enzo's chemistry and heating routines
through specification of photo-heating and photo-ionization rates that
are passed to Enzo's rate solvers.  To that end, the Hydrogen + Helium
photo-heating rate is set to be 
\begin{align}
  \notag
  G &= \frac{c\,E\,\mn_{HI}}{\rhob} 
    \left[\int_{\nu_{HI}}^{\infty} \sigma_{HI}\, \chi_E
    \left(1-\frac{\nu_{HI}}{\nu}\right)\, d\nu\right] \bigg/
    \left[\int_{\nu_{HI}}^{\infty} \chi_E d\nu\right] \\
  \label{eq:photogamma}
    &+ \frac{c\,E\,\mn_{HeI}}{\rhob} 
    \left[\int_{\nu_{HeI}}^{\infty} \sigma_{HeI}\, \chi_E
    \left(1-\frac{\nu_{HeI}}{\nu}\right)\, d\nu\right] \bigg/
    \left[\int_{\nu_{HI}}^{\infty} \chi_E d\nu\right] \\
  \notag
    &+ \frac{c\,E\,\mn_{HeII}}{\rhob} 
    \left[\int_{\nu_{HeII}}^{\infty} \sigma_{HeII}\, \chi_E
    \left(1-\frac{\nu_{HeII}}{\nu}\right)\, d\nu\right] \bigg/
    \left[\int_{\nu_{HI}}^{\infty} \chi_E d\nu\right],
\end{align}
where if {\tt RadHydroChemistry} is set to 1, only the first of these
terms is used.  Similarly, the photo-ionization rates are calculated as
\begin{align}
  \label{eq:kphHI}
  \Gamma_{HI}^{ph} \ &= \ \int_{\nu_{HI}}^{\infty}
    \frac{c\,E_{\nu}\,\sigma_{HI}}{h\nu}\, d\nu 
  \ = \ \frac{c\,E}{h} 
    \left[\int_{\nu_{HI}}^{\infty} \frac{\sigma_{HI}\, \chi_E}{\nu}\,
      d\nu\right] \bigg/ \left[\int_{\nu_{HI}}^{\infty} \chi_E
      d\nu\right] \\
  \label{eq:kphHeI}
  \Gamma_{HeI}^{ph} \ &= \ \int_{\nu_{HeI}}^{\infty}
    \frac{c\,E_{\nu}\,\sigma_{HeI}}{h\nu}\, d\nu 
  \ = \ \frac{c\,E}{h} 
    \left[\int_{\nu_{HeI}}^{\infty} \frac{\sigma_{HeI}\, \chi_E}{\nu}\,
      d\nu\right] \bigg/ \left[\int_{\nu_{HI}}^{\infty} \chi_E
      d\nu\right] \\
  \label{eq:kphHeII}
  \Gamma_{HeII}^{ph} \ &= \ \int_{\nu_{HeII}}^{\infty}
    \frac{c\,E_{\nu}\,\sigma_{HeII}}{h\nu}\, d\nu 
  \ = \ \frac{c\,E}{h} 
    \left[\int_{\nu_{HeII}}^{\infty} \frac{\sigma_{HeII}\, \chi_E}{\nu}\,
      d\nu\right] \bigg/ \left[\int_{\nu_{HI}}^{\infty} \chi_E
      d\nu\right],
\end{align}
where again, the latter two of these rates are not filled unless {\tt
  RadHydroChemistry} equals 3.  Finally, the frequency-integrated
opacity is computed as,
\begin{align}
\notag
  \kappa \ &= \ 
  \left[\int_{\nu_{HI}}^{\infty} \kappa_{\nu}\,E_{\nu}\,d\nu\right] \bigg/
  \left[\int_{\nu_{HI}}^{\infty} E_{\nu}\,d\nu\right] \\
\label{eq:opacity}
  \ &= \ 
  \mn_{HI} \left[\int_{\nu_{HI}}^{\infty}
    \chi_E\,\sigma_{HI}\,d\nu\right] \bigg/
  \left[\int_{\nu_{HI}}^{\infty} \chi_E\,d\nu\right] \\
\notag
  &+ \ 
  \mn_{HeI} \left[\int_{\nu_{HeI}}^{\infty}
    \chi_E\,\sigma_{HeI}\,d\nu\right] \bigg/
  \left[\int_{\nu_{HI}}^{\infty} \chi_E\,d\nu\right] \\
\notag
  &+ \ 
  \mn_{HeII} \left[\int_{\nu_{HeII}}^{\infty}
    \chi_E\,\sigma_{HeII}\,d\nu\right] \bigg/
  \left[\int_{\nu_{HI}}^{\infty} \chi_E\,d\nu\right].
\end{align}
In all of the above, the assumed radiation spectrum $\chi(\nu)$ handle
the change from the original frequency-dependent radiation equation to
the integrated grey radiation equation.




\section{Numerical solution approach}
\label{sec:solution_approach}

We solve these models in an operator-split fashion, using one of two
separate approaches.  In the first approach, {\tt gFLDSplit} performs
chemistry and heating internally; in the second approach {\tt
  gFLDSplit} performs only the radiation evolution while Enzo performs
the remaining physics.


{\bf Approach 1: internal chemistry and heating}

In this approach, we solve the radiation equation
\eqref{eq:radiation_PDE} separately from the more tightly-coupled gas
energy correction and chemistry equations
\eqref{eq:cons_energy_correction} and \eqref{eq:hydrogen_ionization},
which are evolved together.  These solvers are coupled to Enzo's 
existing operator-split solver framework in the following manner:
\begin{itemize}
\item[(i)] Evolve the radiation implicitly in time [{\tt gFLDSplit}].
\item[(ii)] Evolve the gas energy correction and primordial chemistry
  number densities implicitly in time [{\tt gFLDSplit}]. 
\item[(iii)] Project the dark matter particles onto the finite-volume
  mesh to generate a dark-matter density field $\rho_{dm}$ [Enzo];
\item[(iv)] Solve for the gravitational potential $\phi$ using a
  Poisson equation [Enzo];
\item[(v)] Advect the dark matter particles with the Particle-Mesh
  method [Enzo];
\item[(vi)] Evolve the hydrodynamics equations using an up to
  second-order explicit method, and have the velocity $\vb$ advect
  both the Hydrogen I number density $\mn_{HI}$ and the grey radiation
  field $E$ [Enzo]; 
\end{itemize}

The implicit solution approach for step (i) is similar to the one from 
\cite{ReynoldsHayesPaschosNorman2009}; here we describe only enough
to point out the available user parameters, and more fully describe
some additional options available in the solver.




{\bf Approach 2: external chemistry and heating}

In this approach, we solve the radiation equation
\eqref{eq:radiation_PDE} in one phase, and all other Enzo components
are solved in the existing operator-split manner:
\begin{itemize}
\item[(i)] Evolve the radiation implicitly in time [{\tt gFLDSplit}].
\item[(iii)] Project the dark matter particles onto the finite-volume
  mesh to generate a dark-matter density field $\rho_{dm}$ [Enzo];
\item[(iv)] Solve for the gravitational potential $\phi$ using a
  Poisson equation [Enzo];
\item[(v)] Advect the dark matter particles with the Particle-Mesh
  method [Enzo];
\item[(vi)] Evolve the hydrodynamics equations using an up to
  second-order explicit method, and have the velocity $\vb$ advect
  both the Hydrogen I number density $\mn_{HI}$ and the grey radiation
  field $E$ [Enzo]; 
\item[(vii)] Evolve the gas energy and chemistry systems  [Enzo]. 
\end{itemize}

The implicit solution approach for step (i) is identical to the
approach used in \label{sec:solution_approach1}.


In solving the steps (i) and (ii) for both approaches, we first
discretize the equations \eqref{eq:radiation_PDE},
\eqref{eq:cons_energy_correction} and \eqref{eq:hydrogen_ionization}
in space and time before we solve them computationally.  We use a
method of lines approach for the space-time discretization, wherein we
first discretize in space, and then evolve the resulting system of
ODEs in time.  As with the rest of Enzo, we use a finite-volume
spatial discretization, placing all of our unknowns at the center of
each finite-volume cell, and performing all spatial derivatives
through a divergence of face-centered fluxes. 



\subsection{Radiation subsystem}
\label{sec:rad_solve}

We discretize the radiation equation \eqref{eq:radiation_PDE} using a
standard two-level $\theta$-method,
\begin{align}
  \label{eq:radiation_PDE_theta}
  E^n - E^{n-1} &- \theta\dt\left(\nabla\cdot\(D\,\nabla E^n\) - \frac{\dot{a}}{a} E^n -
    c\kappa^n E^n + \eta^n\right) \\ 
  \notag
  & - (1-\theta)\dt\left(\nabla\cdot\(D\,\nabla E^{n-1}\) - \frac{\dot{a}}{a} E^{n-1} -
    c\kappa^{n-1} E^{n-1} + \eta^{n-1}\right) = 0,
\end{align}
where $0\le\theta\le 1$ defines the time-discretization, and where we
have assumed that the advective portion of \eqref{eq:radiation_PDE}
has already been taken care of through Enzo's hydrodynamics solver.
Recommended values of $\theta$ are 1 (backwards Euler) and $\frac12$
(trapezoidal, a.k.a.~Crank-Nicolson).  

Whichever $\theta$ value we use (as long as it is nonzero), the
equation \eqref{eq:radiation_PDE_theta} is linearly-implicit in the
time-evolved radiation energy density $E^n$.  We write this in
predictor-corrector form (for ease of boundary condition
implementation), which we will write as
\begin{align}
\label{eq:linear_system}
  J s = b, \qquad E^n = E^{n-1} + s.
\end{align}
We approximately solve this linear equation for the update $s$,
to a tolerance $\delta$,
\begin{align}
\label{eq:linear_system_approx}
  \| J s - b \|_2 \le \delta,
\end{align}
using using a multigrid-preconditioned conjugate gradient iteration.




\subsection{Gas and Chemistry subsystem}
\label{sec:analytic_solve}

This phase is only performed if {\tt gFLDSplit} is performing the
chemistry and heating solve.  As with the radiation equation, we also
assume that the advective portion of \eqref{eq:hydrogen_ionization} is
taken care of using Enzo's hydrodynamics solvers.  Therefore, since
the remainder of the gas energy and chemistry equations is spatially
local, the solver for step (ii) is performed separately on a
cell-by-cell basis.  To this end, we employ a new implicit-time
version of the {\em quasi-steady-state approximation}.  In this
approach, instead of approximating the solution to the exact ODEs, we
exactly solve approximate ODEs.  Here, if we assume in each equation
that all but the time-evolving variable are held constant throughout
the time step, we could consider the ODE system 
\begin{align}
  \label{eq:energy_correction_qss}
  \partial_t e_c &= -\frac{2\dot{a}}{a}e_c + G(\overline{E},\overline{\mn_{HI}}) - \Lambda(\overline{E},e,\overline{\mn_{HI}}), \\
  \label{eq:hydrogen_qss}
  \partial_t \mn_{HI} &= \alpha^{rec}(\overline{e}) \mn_e \mn_{HII} - \mn_{HI} \Gamma_{HI}^{ph}(\overline{E}),
\end{align}
where $\overline{u}$ denotes a field $u$ that is assumed fixed
throughout a time step.  With this approximation,
\eqref{eq:energy_correction_qss}-\eqref{eq:hydrogen_qss} may be
written as
\begin{align}
  \label{eq:energy_correction_qss2}
  \partial_t e_c &= Q - Pe_c, \\
  \label{eq:hydrogen_qss2}
  \partial_t \mn_{HI} &= a \mn_{HI}^2 + b\mn_{HI} + c,
\end{align}
where $P$, $Q$, $a$, $b$ and $c$ are all constant throughout the time
step.  These may be solved analytically to full accuracy for any time step
size $\dt$.  We write these analytical solvers as
\begin{align}
  \label{eq:energy_correction_qss3}
  e_c(t) &= \text{sol}_e\left(\overline{E},\overline{\mn_{HI}},e_c^{n-1},t\right), \\
  \label{eq:hydrogen_qss3}
  \mn_{HI}(t) &= \text{sol}_{HI}\left(\overline{E},\overline{e},\mn_{HI}^{n-1},t\right),
\end{align}
and couple these together implicitly through defining the nonlinear
equations 
\begin{align}
  \label{eq:energy_correction_iqss}
  f_e(e_c,\mn_{HI}) = e_c^n &- \text{sol}_e\left(\frac{E^{n-1}+E^n}{2},\frac{\mn_{HI}^{n-1}+\mn_{HI}^n}{2},e_c^{n-1},t\right) = 0, \\
  \label{eq:hydrogen_iqss}
  f_{HI}(e_c,\mn_{HI}) = \mn_{HI}^n &- \text{sol}_{HI}\left(\frac{E^{n-1}+E^n}{2},\frac{e^{n-1}+e^n}{2},\mn_{HI}^{n-1},t\right) = 0.
\end{align}

As solutions to these equations can change rather dramatically
(e.g.~ionization states may change by orders of magnitude in a cell in
a single time step), we use a highly robust damped fixed-point
iteration to solve the equations
\eqref{eq:energy_correction_iqss}-\eqref{eq:hydrogen_iqss},
\begin{equation}
  \label{eq:nonlinear_system}
  \begin{split}
  |f_e(e_c,\mn_{HI})| &< 10^{-8}, \\
  |f_{HI}(e_c,\mn_{HI})| &< 10^{-8},
  \end{split}
\end{equation}
at each time step to obtain the updated solution variables $e_c^n$ and
$\mn_{HI}^n$.

We note that if the local thermodynamic equilibrium approximation is
used, then this subsystem only solves for the time-evolved gas energy
correction $e_c^n$.




\subsection{Time-step selection}
\label{sec:dt_selection}

Time steps are chosen adaptively in an attempt to control error in the
calculated solution.  To this end, we first define an heuristic
measure of the time accuracy error in a specific variable $u$ as
\begin{align}
\label{eq:time_error}
  err = \left(\frac1N \sum_{i=1}^N
    \left(\frac{u_i^{n}-u_i^{n-1}}{\omega_i}\right)^p\right)^{1/p}, 
\end{align}
where the weighting vector $\omega$ is given by
\begin{align}
\label{eq:time_weighting}
  \omega_i &= \sqrt{u_i^n u_i^{n-1}} + 10^{-3}, \quad i=1,\ldots,N, \\
  \omega_i &= |e_{c,i} + e_{h,i}| + 10^{-3}, \quad i=1,\ldots,N,
\end{align}
i.e.~we scale the radiation and chemistry change by the geometric mean
of the old and new states, and scale the gas energy change by the new
total gas energy, adding on a floor value of $10^{-3}$ in case any
of the states are too close to zero.  This approach works well when
the internal solution variables are unit-normalized, or at least close
to unit-normalized, since the difference between the old and new
solutions, divided by this weighting factor $\omega$, should give a
reasonable estimate of the number of significant digits that are
correct in the solution. 

With these error estimates \eqref{eq:time_error} for each variable, we
set the new time step size for each subsystem based on the previous
time step size and a user-input tolerance $\tau_{\text{tol}}$ as
\begin{align}
\label{eq:time_estimate}
  \dt^{n} = \frac{\tau_{\text{tol}} \dt^{n-1}}{err}.
\end{align}
Since $E$ and $\{e_c,\mn_{HI}\}$ are evolved separately,
we allow the $\{e_c,\mn_{HI}\}$ solver to subcycle at a faster rate
than the $E$.  We therefore have two time step sizes that we use in
the module,
\begin{align}
\label{eq:FLD_time_estimate}
  \dt_{E}^{n} &= \min\{\dt_{E}^{n},\dt_{CFL}^{n}\}. \\
  \dt_{e,HI}^{n} &= \min\{\dt_{e}^{n},\dt_{HI}^{n},\dt_{E}^{n}\},
\end{align}
where $\dt_{\text{CFL}}$ is the time step size that Enzo's other
routines (e.g.~hydrodynamics) would normally take.  A user may
override these adaptive time step controls with the input parameters 
$\dt_{\text{max}}$ and $\dt_{\text{min}}$. 

We further note that when run in combination with Enzo's hydrodynamics
routines, the {\tt gFLDSplit} solver will limit Enzo's hydrodynamics
solver to have a maximum time step of size 
{\tt RadHydroMaxSubcycles}$\dt_{\text{E}}$.  That said, the radiation
time step will {\em never} be larger than the hydrodynamics time step
size.  As a result, in some physical regimes, the global time step
size will be limited based on the radiation time scale, and in other
regimes it will be limited by the hydrodynamic time scale.




\subsection{Variable rescaling}
\label{sec:variable_rescaling}

In case Enzo's standard unit non-dimensionalization using 
{\tt DensityUnits}, {\tt LengthUnits} and {\tt TimeUnits} is
insufficient to render the resulting solver values $E$, $e_c$ and
$n_{HI}$ to have nearly unit magnitude, the user may input additional
variable scaling factors to be used inside the {\tt gFLDSplit}
module.  Denoting these user-input values as $s_E$, $s_e$ and
$s_{\mn}$, then we may define the rescaled variables
\begin{align}
\label{eq:variable_rescaling}
  \tilde{E} = E / s_E, \qquad \tilde{e}_c = e_c / s_e, \qquad 
  \tilde{\mn}_{HI} = \mn_{HI} / s_{\mn},\qquad 
  \tilde{\mn}_{HeI} = \mn_{HeI} / s_{\mn},\qquad 
  \tilde{\mn}_{HeII} = \mn_{HeII} / s_{\mn},
\end{align}
and the {\tt gFLDSplit} module will use $\tilde{E}$, $\tilde{e}_c$ and
$\tilde{\mn}_{HI}$ in its internal routines instead of Enzo's internal
variables $E$, $e_c$ and $\mn_{HI}$.  If the user does not know
appropriate values for these scaling factors {\em a-priori}, a
generally-applicable rule of thumb is to first run their simulation
for a small number of time steps and investigate Enzo's HDF5 output
files to see the magnitude of the values stored internally by Enzo; if
these are far from unit-magnitude, these scaling factors should be
used. 

In addition, if enabled through the {\tt AutomaticScaling} input
parameter, {\tt gFLDSplit} can automatically update the values of
$s_E$, $s_e$ and $s_{\mn}$.  These are set to be the maximum absolute
value of each internal field (e.g. $\tilde{E}$, etc) from the previous
time step, such that at the beginning of the next time step the
maximum value of $\tilde{E}$ and $\tilde{e}_c$ equal 1, and the
maximum value of the largest of $\tilde{\mn}_{HI}$,
$\tilde{\mn}_{HeI}$ and $\tilde{\mn}_{HeII}$ equals 1.  We note that
this approach \underline{will fail} if the radiation field values
throughout a domain can deplete to 0, since eventually $s_E$ will
decrease in magnitude until $\tilde{E} \approx \frac{0}{0}$.  In these
situations, it is instead recommended to set these scaling factors
manually, and to disable {\tt AutomaticScaling}.




\subsection{Boundary conditions}
\label{sec:boundary_conditions}

As the radiation equation \eqref{eq:radiation_PDE} is parabolic,
boundary conditions must be supplied on the radiation field $E$.  The
{\tt gFLDSplit} module allows three types of boundary conditions to
be placed on the radiation field:
\begin{itemize}
\item[0.] Periodic,
\item[1.] Dirichlet, i.e.~$E(x,t) = g(x), \; x\in\partial\Omega$, and
\item[2.] Neumann, i.e.~$\nabla E(x,t)\cdot n = g(x), \; x\in\partial\Omega$.
\end{itemize}
In most cases, the boundary condition types (and values of $g$) are
problem-dependent.  When adding new problem types, these conditions
should be set near the bottom of the file {\tt gFLDSplit\_Initialize.C}, 
otherwise these will default to either (a) periodic, or (b) will use
$g=0$, depending on the user input boundary condition type.



\section{Concluding remarks}
\label{sec:conclusions}

We wish to remark that the module is not large (one header
file, 15 C++ files, 6 F90 files), and all files begin with the 
{\tt gFLDSplit} prefix.  While we have strived to ensure that the
module is bug-free, there is still work to be done in enabling
additional physics, including Helium/molecular chemistry and more
advanced time-stepping interactions with the rest of Enzo (especially
when ionization sources ``turn on'' abruptly).  

Feedback/suggestions to are welcome.


\bibliography{sources}
\bibliographystyle{siam}
\end{document}
